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ABSTRACT 

We show that for young stars which are still accreting and for which measurements of 
stellar age, t*, disc mass Mdisc and accretion rate M are available, nominal disc age 
idisc = Mdisc/ M is approximately equal to the stellar age at least within the considerable 
observational scatter We then consider theoretical models of proto-stellar discs through 
analytic and numerical models. A variety of viscosity prescriptions including empirical 
power laws, magnetohydrodynamic turbulence and gravitational instability were considered 
within models describing the disc phenomena of dead zones, photoevaporation and planet 
formation. These models are generally poor fits to the observational data, showing values of 
idisc which are too high by factors of 3 - 10. We then ask whether a systematic error in the 
measurement of one of the observational quantities might provide a reasonable explanation 
for this discrepancy. We show that for the observed systems only disc mass shows a system- 
atic dependence on the value of idisc/i* and we note that a systematic underestimate of the 
value of disc mass by a factor of around 3-5, would account for the discrepancy between 
theory and observations. 



1 INTRODUCTION 

Proto-stellar discs form during the star formation process. Much 
of the material falling inwards during the final collapse phase has 
non-zero angular momentum about the newly formed proto-star. 
This angular momentum must be removed in order for the mate- 
rial t o fall further inwa r ds, an d this is the basis of an accretion 
disc jWilliams & Ciez j bOUh ; lArmitag3 ( [2010,) and references 
within). 

The inward falling material flattens into a disc and then vis- 
cous processes, probably driven by a combination of magnetohy- 
drodynamic turbulence and gravitational instabilities, redistribute 
the angular m omentum and allow accretion on to the central star 
( |PringIelll98ll) . Within a few tens of Myr the gas in the disc has ei- 
ther accreted on to the star or been blown away by stellar radiation, 
while dust grains have aggregated, leaving behind a debris disc and 
potentially a planetary system. 

Having a detailed model of these systems is vital to under- 
standing both the late stages of star formation and planet forma- 
tion processes. At present discs are poorly understood and many 
of their parameters are only weakly constrained. These loose con- 
straints stem from the difficulties involved in observing systems 
which subtend just ~ 1" on the sky and for which detection of 
their main constitue nt, molecular hydrogen, is not currently fea- 
sible ^Andrews et al. 2009). Despite this, such discs are often ob- 
served around T Tauri stars through the IR excess they produce 
and emission from their dust in the millimetre and sub-millimetre. 
As the temperature of the disc changes in the radial direction away 



from the central star (as does the disc structure), the spectral energy 
distribution can be used to p robe the disc structur e without it need- 
ing to be spatially resolved jSeckwith et alj[l990D . However, large 
uncertainties remain in many of the disc properties, principally due 
to the assumptions required in order to proceed in this manner. 

In general the equations describing the evolution of accretion 
discs are non-linear and require numerical treatments, but if simpli- 
fying assumptions about the fo rm of the viscosity are made, then 
analytic solutions are possib le ( lvonWeizsackerlll948l : iLiistll 19521: 
iLvnden- Bell&Pringl3ll974h . These solutions give the total disc 
mass decaying as Mdisc oc t~'' , where a > 0. Differentiating 
this leads to the accretion rate onto the central star behaving as 
M oc t-^'^+^K For this reason, as well as on purely dimensional 
grounds, one can define an age indicator, namely the disc age fdisc 



^disc 



Mdis, 
M 



(1) 



Disc age is the characteristic accretion timescale of the disc at any 
given time. For simple viscous disc models disc age is proportional 
to the age of the accretion disc with a proportionality constant of 
order unity. Therefore we regard disc age as a proxy for the age 
of the system, which in principle should be related to the age of 
the central star, t,. Disc age is measurable as both Mdisc & M 
are observable quantities, and can be measured from independent 
observations (see §2). 

In Section 2 we summarise the techniques involved in mak- 
ing the observational measurements and display the observational 
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Figure 1. Observations of young stellar objects with accretion discs. Nom- 
inal disc age, tdiso (y-axis) plotted against stellar age, tt, (x-axis) on a 
logarithmic scale. The line tjisc = is included for comparison. This 
line, tdisc = t*, is a reasonable fit to the data, although there is a high de- 
gree of scatter, consistent with observational uncertainties. The horizontal 
error bar is a factor of 3 (the error estimate of stellar ages). The vertical eiTor 
bar is a factor of 5 from the combination of the error estimates of accretion 
rate and disc mass (factors of 3 and 4 respectively). The detection limits on 
accretion rate and disc mass are approximately 10"^'' Mqyt-'^ & 0.001 
Mq respectively. This may lead to preferential exclusion of data points in 
the upper right of the plot, as disc mass generally decreases with age. The 
minimum stellar age cut off is ~ 10^ Myr 



data in the — tdisc plane. In Section 3 we then investigate vari- 
ous simple theoretical models of accretion disc evolution and show 
that these models tend to produce values of fdisc that are too high by 
factors of 3 - 10. We show that taking account of various complica- 
tions such as photoevaporation and planet formation do not account 
for this discrepancy. In Section 4 we discuss these results and in- 
vestigate the possibility that the discrepancy might be accounted for 
by a systematic error in one of the observed quantities, concluding 
that a possible solution is that the disc masses are systematically 
underestimated by a factor in the range 3-5. 



by the presence of a circum-stellar disc. Accretion onto the star cre- 
ates a blue excess and emission from the disc creates a red excess, 
hence making the spectrum both too blue and too red ( I Hillenbrand 
[2002). In addition starlight may be either reflected off, or extin- 
guished by, a disc, either increasing or decreasing the observed lu- 
minosity. 

Also, there is uncertainty about which theoretical tracks are 
the most relevant, and different authors use different tracks, often 
leading to subs tantially different ag e estimates, particularly for the 
youngest stars. iHillenbrand |2009|) estimates the errors to be ap- 
proximately a factor of 3 in either direction. The fractional errors 
are likely to reduce as the stars become older. 



2.2 Stellar Accretion Rate M 

If the distance and spectral type are known, stellar radii can be 
found from luminosities using the Stephan-Boltzmann law. Mate- 
rial is then assumed to fall from the inner edge of the disc (usually 
taken as Rin ~ 5i?, ) onto the surface of the star, releasing its grav- 
itational potential energy as accretion luminosity, given by 



Lac 



(2) 



where Af, is the stellar mass derived from evolutionary tracks, M 
is the accretion rate, R, is the stellar radius. 

The characteristic temperature of this emission is of the order 
of 10,000 K, which is greater than the stellar emission tempera- 
ture at a few thousand K. The UV excess is measured and taken to 
be entirely due to the accreting material, giving the accretion rate 
through the above relation. 

Uncertainty in the distance to the sources could be a substan- 
tial source of systematic error. Converting from stellar luminosity 
to radius introduces one factor of distance, ci,and converting from 
observed accretion flux to luminosity adds two more, meaning the 
final value depends on d'^. The other key source of error is the stel- 
lar mass derived from evolutionary tracks. Final errors in the ac- 
cretion rate were taken to be a f actor of 3 ( iHartmann et al.|[l998l : 
ICalvet. Hartmann & Strori]|2000h . 



2 OBSERVATIONAL DATA 

In principle the three properties required to compare the disc age 
fdisc with the stellar age for a given system (accretion rate, disc 
mass, and stellar age, taken as the system age) are all measurable. 
The requisite observational data, taken from a variety of sources, 
are presented in Figure[T]and Table[T] In this section we give a brief 
outline of the procedures behind the measurement of each property. 
For a more detailed description we refer the reader to the referenced 
sources. 

2.1 Stellar Age 

Spectroscopy is used to determine the spectral type and effective 
temperature of each star. Combined with the measured luminosity 
this allows the star to be located in the Hertzsprung-Russell dia- 
gram (hereinafter H-R diagram). From here theoretical pre-main- 
sequence evolutionary tracks are used to convert the H-R diagram 
position to a stellar age and mass. 

As with all astrophysical sources the light from these systems 
is reddened by intervening material and this must be corrected for. 
However, the extinction correction process may well be confused 



2.3 Disc Mass Mdisc 

Most of the mass of a proto-stellar disc is at large radii and is 
contained in molecular hydrogen. However this is a symmetric 
molecule with no permanent electric dipole, making it difficult to 
detect. Thus the discs are usually detected through IR emission 
from warm dust in the inner disc, or through mm or sub-mm emis- 
sion from cool dust in the outer disc. 

To estimate the disc mass from the IR fluxes it is neces- 
sary to have a simple model of the disc structure in mind. It is 
usual to assume surface density structure and a temperature dis- 
tribution of the forms T.{R) = EoR'^ and T{R) = Toi?"' re- 
spectively. Then by assuming a form for the dust op acity (usually 
= 0.1(i^/10^^ Hz) g cm"2 ( lBeckwithetalJl99(]|) ) and that the 
dust emits as a black body at each radius, the free parameters p, q. 
To and Eq can be adjusted to achieve the best fit to the observed 
spectrum. Since the IR flux comes from warm dust, this gives the 
total dust mass in the disc. The total disc mass is assumed to be 100 
times this value. 

The primary source of error here is thought to be the assumed 
dust opacity. Grains are expected to grow and alter the form of the 
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Table 1. Observational data on young stellar objects with accretion discs 
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Column 1: Object name; Column 2: logj^Q of stellar age estimate (yr); Column 3: Disc mass 
estimate (Mq), the dust masses from source R were averaged between q = 2.5&;3, and multiplied 
by 100 to give total disc mass; Column 4: log^o of stellar accretion rat e (Mq yr~^); Co lumn 
5: login of disc age (= MHj .sr./A-f) in years; Column 6: S ource s R = iRicci et al ] fcOlOl) . A7 
= lAndrews^ & Williamsj jioOTI) A9 = [An drews et al.l b009|) H = iHartmann et al.' )l998| ) V = 
Various jw hite & G he3 l200ll; iNaiitaTs trom & MuzeroM l2007t [Andrews & WilliamT l2003; 
iMuzerolleet al.l2003tlG^llbring et all 19981) . 



opacity (as well as the gas to dust ratio), but this process is not well 
understood and a simple form must be used at present. 

Taking the ISM average dust to gas ratio is probably a reason- 
able value for the initial disc, but this is likely to change as material 
aggregates, maybe forms planets and gas is preferentially lost in 
winds. However, it is thought that the lack of treatment of grain 
growth will decrease the mass estimates, as larger bodies are ne- 



glected, where as using the ISM dust to gas ratio will over estimate 
the mass, so there should be some degree of cancellation between 
these two sources of error. The overall errors i n total disc mass were 
typic ally estimated to be a factor of 3 or 4 ( [Andrews & Williams! 
i200g[ ; [Ricci etan[201(]|) . 
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2.4 Comparison of disc age and stellar age 

As mentioned above we present in Figure [T] the observational es- 
timated values for both fdisc and t,, for those stars for which we 
were able to discover estimates of mass, accretion rate and age in 
the literature. We also indicate the sizes of the estimated errors in 
both these quantities. In Figure [T] we also plot the line t^isc ~ t*- 
Given the size of the estimated errors this line seems to provide a 
reasonable fit to the data, although there is considerable scatter. 



3 THEORETICAL MODELS & RESULTS 

In this Section we investigate the extent to which theoretical models 
of accretion disc evolution agree with the observational estimates 
in the i, — tdisc plane. ^ 

For an axisymmetric matter distribution with surface density 
E(i?) and viscosity viR) in Kelperian orbits ar ound a central point 
mass, the surface density evolves according to l iPringlelll98ll) : 

9E _ 3__9_ 

dt ^ RdR 

In general v is likely to be a function of local disc variables such as 
E and R, making the equation non-linear and requiring solution by 
numerical means. There are however special cases where analytic 
results can be found. 



R 



1/2 d 
dR 



(3) 



3.1 Similarity Solutions 

If u is taken to be a power law function of radius R alone then 
the equation for E is linear and an analytic solution to the initial 
value problem is possibl e by means of a Green' s function ( LUst 
ll952l;lLvnderi^Bell & Pringldll974|Pringlelll98ll) . For viscosities 
only dependent on a power of the radius (v oc 7?") the solutions 
which correspond to zero torque at the origin, so that the total an- 
gular momentum of the disc is conserved, but accretion takes place 
at small radius, take on similarity form at late times: 



C 



exp 



(4) 



where r = R/Ro, t = t/{to + 1), p is a simple algebraic expres- 
sion in n, and C, -Ro and to are constants. 
Thus at late times we find 



Mdisc oc t 
where 

1 



2(2 



(5) 



(6) 



The requirement that a > implies that n < 2. 

In the more general case in which u oc E""/?" with m > 
— 1 for viscous stability (Pringle 1981), there are also similarity 
solutions. They are of of the form (Pringle 1991) 



(7) 



where Eq, Ro and to are model dependent constants, and a, /?, 7, 
7], fi and k are simple algebraic expressions in m and n. 
At late times these solutions also give 

^ For clarity t will be assumed equal to henceforth. Stiictly they are not 
equal as t is the time since the beginning of the simulation, not the birth of 
the star. This approximation improves with age. 



Table 2. Opacity regimes (n = l — 12) for the opacity law approximated 
as a series of power law segments k = KopT'' cm^ Maximum 
tempe rature Tmax is g iven for regimes independent of density i B ell & Lii] 
1994; lBell et alj|l997l) . The value of u relevant to the mass decline in the 
corresponding similarity solutions (Equations[8]&|9) is also given for those 
regimes n which correspond to a thermally and viscously stable disc. 
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Mdisc oc t 
where now 



1 



2n 



(8) 



(9) 



^disc 



(10) 



5m -I- 4 

These similarity solutions correspond to all the mass initially be- 
ing at the origin at time t = Q, but with a finite amount of angu- 
lar momentum. Mass is then accreted at the origin, with the disc 
spreading to conserve angular momentum. Eventually all the mass 
accretes onto the central object and all the angular momentum is 
carried away to infinity by an infinitesimal amount of the mass. 

Because for these solutions Afdisc oc t^" , which can be dif- 
ferentiated to give M oc at~^-^^'^^ the disc age for these solutions 
is given by 

A'/disc _ t_ 

This represents a straight line parallel to tdisc = t* in the log(t*) - 
log(idiac) plane, but translated vertically by a distance log(l/a). 



3.2 Basic Models 

We first investigate the properties of some simple models for disc 
evolution. 

We solve the equations governing the evolution of surface den- 
sity using an explicit first order finite difference method on a radial 
grid of 100 points, logarithmically spaced in radius R. We take the 
inner disc boundary radius to be i?in = 5 Rq and the outer bound- 
ary radius i?out = 1000 AU. The outer boundary is chosen so that 
is has little or no effect, except at very late times log(t/yr) > 7.5, 
when in some models the disc touches the outer radius. The bound- 
ary conditions were zero torque (E(i?in) = 0) at the inner bound- 
ary, and zero radial velocity at the outer boundary. We take the 
stellar mass to be constant throughout (Af, = 0.75 Mq), and so 
ignore the mass of any added disc material. Typically we assume 
matter is added to the disc at around Ro — 10 AU with the local 
disc temperature and specific angular momentum. Initially all the 
disc mass is contained in a narrow Gaussian of standard deviation 
1 AU centred on Rq, with an initial total disc mass of 0.2 A/q. 
These conditions hold for all of the following simulations unless 
explicitly stated otherwise. 
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3.2.1 Power-law viscosity 

We first consider a simple viscosity of the form 

u = uo(R/Ro), 



(11) 



witfi uq and Rq constants. Tliis viscosity is implied by the common 
assumption that the disc surface density is of the form E oc R~^. 
The loci of the solutions in the — tdisc plane are given in Figure|2] 
where we have also illustrated the effect of varying the initial disc 
radius Ro, for given uq. This is, of course, equivalent to varying 
the size of the viscosity z^o for fixed initial disc radius. At late times 
the solutions approach the similarity solutions, and so idisc/t* ^ 
1/a — 2, since here n = 1. Note, however, that the approach 
to these solutions is from above, since at early times the accretion 
rate is lower for a given disc age, since the disc takes time to fully 
extend to small radii. 



3.2.2 a viscosity 

A more physically motivated approach to the viscosity is to assume 

(12) 



originally from lShakura & SunvaevI ( Il973h . where Cg — TZT/ ^ is 
the sound speed squared, Q{R) is the orbital angular velocity, TZ 
is the gas constant and fi is the mean molecular weight, here taken 
to be 2.3 (appropriate for a molecular disc). The quantity a is the 
standard dimensionless measure of the viscosity. 

In order to implement the a prescription a local model of the 
disc structure is required. We take the usual one-zone model, which 
assumes thermal and hydrostatic equilibrium in the local disc an- 
nulus. The heating rate (per unit area, for each side of the disc) in 
the local annulus due to viscous dissipation is 



Q+{R) = lu{R)E{R)n{Rf 

O 



(13) 



and this is set equal to the local cooling rate due to blackbody emis- 
sion 



=C^SBTe(i?)^ 



(14) 



where here ctsb is the Stefan-Boltzmann constant and Tc is the 
effective surface temperature of the disc. 

To specify the surface temperature To, assumptions regarding 
the vertical structure must be made. We adopt the simplistic ver - 
tically averaged approach from lArmitage, Livio & Pringld j200lh . 



givmg 

3r 



(15) 



where Tc is the temperature on the central disc plane. 

This, together with hydrostatic equilibrium perpendicular to 
the disc plane, which gives the disc scale-height 



H 



(16) 



and the assumption that disc density (on the central plane) is related 
to surface density by pc = E/2J/, form a set of equations which 
specify the disc structure. We note that then the optical depth r is 
given by 



T = -Y;k{pc,Tc) 



(17) 



Thermal equilibrium is a reasonable assumption if the ther- 
mal timescale is much smaller than the viscous timescale 



jPringle. Verbunt & Wadel 19861) . Preliminary tests with a more de- 
tailed thermal energy equation showed this approximation was rea- 
sonable, and so it was adopted as it was less computationally de- 
manding. 



3.2.3 Single Opacity 

We first computed disc evolution using a single opacity which is 
likely to be representative of those regions of the disc which con- 
trol the ev olution, taking k = 0.02r'' ' ^ cm^ g~^, independent of 
density p ( iBell & Lin|[T994l : iBell et alj|l997h . The initial disc was 
taken to be an annulus at radius _Ro = 10 AU, as described above. 
The solutions in the — tdisc plane are shown in Figure [3] for 
three different values of a = 0.1,0.01,0.001. A mean value of 
a — 0.01 is known to give d isc evolutionary times cales in line 
with known properties of discs jHartmann et alj 19981) . 

As can be seen from the figure, at late times the discs all tend 
to the same similarity solution, which in this case gives tdisc/i* — > 
1/(7 = 4.44. As before, the loci all approach the eventual power 
law locus from above, because at early times the disc has yet to 
extend to the centre, and so the accretion rate is lower than that 
predicted by the similarity solution. The approach occurs later for 
smaller values of a. 



3.2.4 Full Opacity 

We next computed the evolution of such a disc using a set of 
opacities given in Table 2. Also given in Table 2 are the values 
of 1/(T corresponding to each temperature range (those opacities 
which give rise to discs which are thermally/viscously unstable are 
marked with an asterisk in this column.). In Figure |4]we compare 
the loci of a disc using the full opacity table with the disc with 
only a single opacity. We have chosen one of the larger values of a, 
0.01, to illustrate the difference in the similarity solution part of the 
loci. Even for the full opacity, the locus tends to a constant value 
of tdisc/t* ~ 10, because at late times the disc evolution is being 
controlled by a single low-temperature regime in the disc at large 
radius. 



3.2.5 Full Opacity with Delayed Infall 

As we have seen, although the tdisc = t* line provides an adequate 
fit to the data, the standard theoretical disc models with realistic 
physical opacities tend to produce loci in the — tdisc plane which 
lie factors of 4 - 10 above it . Altering the evolutionary timescale 
of the discs, by changing the initial disc radii or the size of the 
viscosity merely exacerbates the problem. 

Rather than start with all the mass in the disc at time t = 
we have therefore experimented with the more realistic assumption 
that the mass is added to the disc on some finite timescale. Usual 
estimates for the timescale of disc formatio n are in the region of 
10^ years ( iLarsoiJll969i ; lLin & Pringldl 19901) . 

In Figure |5] we show two additional examples of the result- 
ing behaviour For the first, we start at t = with a disc annu- 
lus of mass AI = O.O2A^0, and then allow an additional infall of 
O.OIA/q onto the annulus to start at a time of t = 10^' yr. The ini- 
tial accretion rate at that stage is taken to be 1 x lO'^^M© yr~^ and 
decays exponentially with a timescale of 10® yr For the second, the 
initial disc annulus of O.2Af0 is fed at an initial rate of 3 x 10® 
yr~^, which falls off exponentially on a timescale of 10® yr 
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Figure 2. Increasing the radius of the centre of the initial mass distribution 
has a knock on effect on the viscous timescale. With v cc R the viscous 
timescale is proportional to R, so as the initial disc radius Rq is doubled in 
each consecutive model they should be initially separated by approximately 
log 2 along the horizontal axis, but still tend towards the same similarity 
solution at late times. 
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Figure 4. Introducing a more physical opacity has little overall effect other 
than to change the similarity solution the model tends to at late times. This 
raising of the disc age means that this, the most physical of the basic models, 
cannot explain the majority of the observation, which largely fall below the 
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Figure 3. Varying the alpha parameter changes the strength of the viscosity, 
altering the viscous timescale, and hence the time taken to converge to the 
appropriate similarity solution, which is the same regardless of the value of 
alpha. 

As can be seen from Figure [S] late or delayed infall does in- 
deed have the effect of rejuvenating the disc, and so bringing the 
loci in the t* — tdisc plane below the locus of the similarity solu- 
tion. Also shown in Figure|5]are the observational data points. It is 
evident that even late or delayed infall does not really have much 
impact in trying to bring theoretical disc predictions into line with 
the data. 



3.3 Complications to the basic model 

In the above we have considered a basic model for disc evolution 
in which the matter is added to the disc in a simple fashion and the 
disc is allowed to evolve using the standard a-prescription for the 
viscosity. However, in reality disc evolution is more complicated 
that this. We consider here three additional complications. First we 
consider the model of a layered disc, in which the standard vis- 
cosity mechanism is rendered inoperative if the disc is too cool, 
so that the magneto-rotational instability cannot give rise to MHD- 
turbulence. Second, we consider what happens in the later stages of 
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Figure 5. Comparison of three alpha prescription discs (with full opacity), 
the first has no mass infall, the second is an initially low mass (0.02 Mq) 
disc with 0.01 Mq of mass infall around t = 1 X 10** yr, and the third be- 
gins with mass infall at 3 X 10"'' Mq yr~^ with a fall off time of 10^ yr. 
The mass infall has a rejuvenation effect on the disc age, but even doubling 
the mass of a disc in a small fraction of its lifetime does not produce a large 
enough effect to to drop below the line (disc = Also plotted are the 
observational data. The majority of the observations are out of reach with 
these models alone, as under reasonable constraints they cannot be made to 
drop below the line (disc = 



the disc's evolution when the disc is evaporated by radiation from 
the central star. And, third, we consider the formation of a planet 
which is massive enough to affect disc evolution. We treat each of 
these effects in fairly idealised form in order to illustrate the likely 
magnitude and direction the effect has on the loci in the — t^isc 
plane. We consider each in turn. 



3.3. 1 Layered Disc - Model 

Proto-stellar discs are thought to have episodic accretion histories 
tenoch et al1l2009h . Such accretion outbursts have been observed 
and are known as FU Orionis events jHartmann & KenvoiJ[T99^ . 
One possible model for this type of phenomena is that of a layered 
disc. 
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The model used here is primarily based on that of 
lArmitage. Livio & Pringld ( 1200 ih . Magnetohydrodynamic (MHD) 
turbulence is assumed to be the process driving viscosity (at small 
radii), but to be operative this requires the gas to be sufficiently 
ionised in order for it to interact with magnetic fields. If the mid- 
plane temperature Tc, in the disc is above a critical value Tcrit = 
800 K then the gas is assumed to be MHD 'active'. The top and 
bottom surface layers (of surface density Siaycr ~ 100 g cm~^) of 
the disc are always assumed to be active as they are kept ionised by 
cosmic rays and/or stellar radiation. Thus an active surface density 
is defined by 



E for Tc > Tcrit 
2Eiaycr for Tc < T 

c 



(18) 
(19) 



Assuming the disc is Keplerian (Q = ^JGM, / R^), the modified 
surface density evolution equation becomes 



9E _ 3__9_ 

dt ^ RdR 



(20) 



where E represents increase in surface density from infalling ma- 
terial. 

Thermal equilibrium and the simplistic vertical structure used 
in i ]3.2.2l are assumed as before. Finally the complete description of 
the viscosity is two-fold. The regions active to MHD turbulence are 
described using the Q-prescription described in i]3.2.2| but in addi- 
tion the model also takes into account the effects of self-gr avity us- 
ing a prescription dependent on the Toomre Q parameter dToomrd 
Il964h . 



Q 



ttGE' 



(21) 



When Q drops below some critical value, taken here to be Qcrit = 
2, the effects of gravitational instability are deemed to be im- 
portant and are accounted for by an addit ional effective Qgrav 
parameter in the Q -viscosity prescription jLin&Pringlll 19871 : 
[Armitage. Livio & P ringle 2001). 



= 0.01 



Qcrit 



- 1 



for Q < Qc 



= otherwise. 



(22) 



(23) 



This complete model results in a disc structure where material at 
small radii (< 0.1 AU) is all highly ionised and MHD active, be- 
yond this region out to a few AU only the surface layers are MHD 
active and there is a 'dead zone' in the disc mid-plane where mate- 
rial builds up. Finally, in the cool outer part of the disc gravitational 
instability is the principle mechanism for angular momentum trans- 
port. 



3.3.2 Layered Disc - Results 

Introducing a layered disc structure has a dramatic effect on the 
evolution of the disc. Figure |6] compares layered and non-layered 
Q-prescription models, demonstrating that the forms of disc ages 
have almost no resemblance to each other. The infall models have 
material added to the outer disc (at 10 AU) at an initial rate of 
3 X lO"'' Mq yr~\ with a fall off time of 10^ yr. The layered 
disc with infalling material initially has longer and more frequent 
accretion outbursts as the dead zone is replenished with mass more 
rapidly than in the model without infall. The difference becomes 
less pronounced at late times (when the infall has decayed away), 
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Figure 6. Comparison of unlayered and layered discs with and without in- 
falling mass. The layered discs are clearly different to the unlayered disc. 
They spend most of their lifetime at high disc age, with intermittent accre- 
tion outbursts dropping them briefly to low disc age. In this way they cover 
almost the complete range of observed disc ages, but are unlikely to be able 
to return the same distribution. 



and is mainly noticeable through the higher disc mass, resulting in 
higher disc age. 

This demonstrates that infalling mass actually has an aging ef- 
fect on the disc, that is, it makes its disc age greater. This is some- 
what counter-intuitive as one would expect higher mass discs to be 
younger and less evolved. The reason for this difference is that un- 
like in unlayered models the total mass of the disc does not play a 
significant role in determining the stellar accretion rate during the 
quiescent periods. The layered structure means much of the extra 
mass is held up in the dead zone and does not contribute to the ac- 
cretion rate (at that time). Thus, the disc age increases because the 
mass increases but the accretion rate does not. 

Outbursts are triggered when the dead zone becomes so dense 
that it is gravitationally unstable {Q < 2). It then becomes viscous 
due to gravitational instability, which releases energy, quickly heat- 
ing the material in the dead zone above the critical temperature for 
it to be MHD active and then the whole zone rapidly accretes onto 
the star. These outbursts drastically alter the disc age of the system 
and are the only model which enters the region below the disc age 
line. 

The loci of these models in the i, — tdisc plane cover almost 
the entire range of the observed disc ages, but because they spend 
most of the time in quiescence, i.e. at low accretion rates, in a 
time-averaged sense the loci are too high to be compatible with 
the data. In order to match the data there would need to be a very 
significant over-detection of outbursting systems, because most of 
the systems are observed at low disc ages, not at high disc ages, 
where a layered disc would spend most of its lifetime. The systems 
undergoing outbursts may not be in that state for a long period, 
but they will be far brighter in the UV and hence easier to detect. 
Thus, it may be the case that the sample displayed is not represen- 
tative of the true distribution. If this were true one would expect 
the observed distribution of accretion rates to be bimodal, with a 
significant number of the systems having very high accretion rates 
(> 10"'' Mq yr~^) and the rest grouped around 10"** Mq yr~^. 
This bimodality should also be reflected in the accretion rates of 
systems above and below the line tdisc ~ t* if over detection of 
outbursts were to explain why so many systems have low disc ages. 
However, inspection of the data indicates that this is not the case. 

The layered disc model does however offer a possible expla- 
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nation for systems observed far above the disc age line. Many of 
the disc ages higher than 10^ years in the observational data are 
that high, in part, due to the low accretion rates in those systems. 
A layered disc spends most of its lifetime in this state, where mass 
is caught up in a dead zone and the system has a low quiescent 
accretion rate, resulting in a high disc age. 



3.3.3 Photoevaporation - Model 

Accretion discs around young stars must undergo some form of dis- 
persal as they are only observed to persist for the order of 10 Myr. 
A favoured model for this dispersal is photoevaporation, where ion- 
ising radiation from the central star causes material in the outer disc 
to be blown off in a wind. The model used here is based on that of 
IClar ke. Gendrin & Sotomavor (2001). 

Material on the surface of the disc is heated by FUV radia- 
tion from the central star. This raises the temperature of the outer 
disc atmosphere and beyond a certain radius can give it sufficient 
ener gy to escape. This radius is know as the gravitational radius 
( iHol lenbach e t al..,1994) . and is given by 



(24) 



where Csi is the sound speed in the photoionised gas. Outside of 
this radius the surface density loss rate due to photoevaporation is 



= 2csino(ii)mH, 



(25) 



where mn is the mass of hydrogen and no [R) is the number den- 
sity at the base of the outflow, which is parameterised by 



Ft 

no{r) = no(i?g) ( 



where 



5/2 



(26) 



Mi^.)^ 5.7X10^ 

and $ is the flux of ionising photons emitted by the central star. 

Although there ar e now more sophisticate d phot oevap- 
oration models (e.g. lOwen. Ercolano & Clarkd feOllh and 
[Alexander & Armitagel ( l2007h ) and there is debate about whether 
FUV, EUV or X-rays play the most important role, all the models 
exhibit essentially the same phenomenon, with a gap opening in 
the disc at late times and the inner disc being rapidly accreted. 
As here we are principally concerned with ascertaining the effect 
of such phenomenon on the disc age of a system this model is 
sufficient for our current purposes. 



-3/2 



127) 



3.3.4 Photoevaporation - Results 

Photoevaporation from the disc has little effect for most of its life- 
time, but at late times when the evaporation rate becomes compa- 
rable to the accretion rate through the disc, it causes a dramatic 
change. A gap in the disc opens up around the gravitational radius 
(Equation I24t. this prevents the inner disc from being replenished 
by material from the outer disc, where most of the mass resides. 
The inner disc is rapidly accreted, after which accretion stops and 
the outer disc is left to gradually evaporate. Thus at late times the 
accretion rate drops to zero and the disc age tends to infinity (Figure 
IT}. Therefore, although this may be an important process allowing 
disc dispersal, it does not help to explain any spread in the disc age. 
Photoevaporation has almost no detectable effect on the disc age 



No Evap 
Evap 

Disc Age — Ago 



4.5 



5.5 



6.5 



7.5 



log(«./yr) 



Figure 7. The effect of photoevaporation on the disc age of a power law 
viscosity disc, with initial disc mass of 0.02 A/q. The only noticeable ef- 
fect occurs rapidly at late times when the photoevaporation opens a gap in 
the disc and the inner disc rapidly accretes, resulting in the accretion rate 
dropping to zero while there is still a finite mass in the outer disc, hence the 
disc age tends to infinity. 



for the majority of the disc's lifetime, and the period where the disc 
age increases significantly is brief. 

This is a simplified model and the form of the evaporation 
would in fact change once the inner disc is removed, but by this 
time the disc age will already be large, and hence consideration of 
a more complex photoevaporation model is unlikely to illuminate 
further effects on the disc age. Howev er, these 'old' discs m ay still 
undergo small outbursts of accretion jMurphvet i5]|201ll) . on the 
timescale of months down to days. It is possible that observations 
of the disc age of such systems could change significantly from 
survey to survey, but this would of only explain scatter on the upper 
extreme of the stellar age range observed, leaving the majority of 
the observations still in need of an explanation. 



3.3.5 Planet Formation & Migration 

It is widely believed planets form in the discs around young stars. 
Furthermore, since many Jupiter mass planets have been observed 
at small radii where it would not have been possible for them 
to form, massive planets are also expected to migrate inwards 
from their formation radius by interaction with the disc. To model 
this, w e use a simple formulation (based on^ Alexander & Armitagel 
( l2009f) ). which considers the effect of planets of around a Jupiter 
mass on the material in the disc, as the planet accretes material and 
migrates. 

Massive planets exert a torque on material in their vicinity and 
thus alter the surface density evolution equation. The evol ution due 
to both planet and disc is described by jLin & PapaloizouH l986l) 



9E 
dt 



Rdr 



3R 



1/2 



d_ 

OR 



{uT.R 



l/2\ 



2AEi?^/2 



,(28) 



where Ecvap is the photoevaporation as described in the previous 
section and A is the torque exerted by the planet, given by 



A = - 



g^GM. / R 



2R 



-7— for i? < a 



A : 



q^GAU ( a 



2R 



A, 



for R > a, 



(29) 



(30) 



Accretion Disc Age and Stellar Age 9 



4 - 



No Planet 

tform = 1-5 X lO" yr 

Disc Age — Ago 



4.5 



5.5 



6.5 



7.5 



log(«./yr) 



Figure 8. Formation of a 0.5 Mjup planet at t = 1.5 X 10 yr in a 

power law viscosity disc with photoevaporation. The initial disc mass was 
0.03 Mq to ensure the planet had a noticeable effect on the disc but did 
not simply fall into the star. The planet effectively 'pins back' the material 
at larger radii than its orbit, this slows the accretion rate and hence raises 
the disc age. 



where q = Mp/A/» (Mp is the planet mass), a is the semi-major 
axis of the planet's orbit and Ap = max(Ji', jr — a\), where H is 
the scale height of the disc. 

In this process angular momentum is transferred between the 
planet and the material in the disc, this causes the planet to migrate 
at a rate given by 



da 



47r 



— / RAT,dR. 



(31) 



The process of angular momentum exchange also acts as an effec- 
tive repulsive force away from the planet because angular momen- 
tum is removed from material inside the planet's orbit and added 
to material outside. Thus, the planet opens up a gap in the disc. 
In this model material is allowed to flow across the gap, some of 
which the planet will accrete. The efficiency of this process, e, is 
approximated by the fitting function 

1/3 



6(Mp) 



= 1.67 



exp 



+ 0.04, 



(32) 



Mjup J yi.5 Mjup ^ 

where £max is taken to be 0.5 I Veras & Armitagj 20041: 
D' Angelo. Hemiing & KlevI l2002l : iLubow. Seibert & Artvmowic j 
19991) . Giving the accretion rate across the gap as 



-A'/p 



1 + e 

and the accretion rate on to the planet as 

Mr, = eA'fdisc, 



(33) 



(34) 



where Afdisc is the accretion rate through the disc well outside the 
planets orbit. This is calculated by assuming a steady state model 
at 3 times the planet's orbital radius, giving 

Mdisc = 37n/(3a)E(3a). (35) 



3.3.6 Planet Formation & Migration - Results 

When a massive planet forms it holds back material in the outer disc 
and prevents it from flowing to the inner disc as rapidly as it other- 
wise would. This in turn lowers the accretion rate onto the star, and 
thus increases the disc age. After it has formed, the planet accretes 
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Figure 10. Compares the most physical basic models from i|3.2l with where 
the observations would he if all the disc masses were increased by a factor 
of 5. This raises many of the observations into the region above the line 
^disc = i*. where they can potentially be explained by one of the models 
studied, within the observational scatter. 



mass and migrates inwards. In the case we compute(Figure[8](, the 
final planet mass was Mp ~ 0.78 Af.jup (initially 0.5 Mjup) and 
it migrated inwards ending up at ~ 5.2 AU from 5.6 AU. After 
a period of rapid change immediately post formation the disc age 
settles down to another, higher disc age, similarity solution which 
is eventually truncated by photoevaporation (described in previous 
subsection). This illustrates the basic effect the formation of a mas- 
sive planet on the disc age. But the effect of the planet is always 
be to increase the disc age by slowing the accretion rate. Thus the 
effect of planet formation is to push the disc age of systems to lie 
further above the line tdisc ~ t*- 



4 DISCUSSION 

We have considered the predicted evolution of theoretical accretion 
discs in the i, — tdisc plane, where 

Afdisc 



M 



(36) 



and have compared the results with observational data, using the 
estimated age of the central star, t,, as a proxy for the age of the 
disc, t. We noted that, despite large scatter, and large uncertainties 
in the observational estimates, the observational data are approxi- 
mately fit by the relationship fdisc ~ t* (Figure[Tll with a scatter of 
around an order of magnitude above and below this line. 

We showed that if the discs are modelled using viscosity pre- 
scriptions which are either a simple power law or the standard alpha 
model, then the loci in the — tdmc plane quickly evolve to a line 
of the form tdisc/t , = 1 /a. For the powe r law of the form p cc R, 
we find 1/(T = 2 jHartmann et alj[l998l) . but for more physically 
realistic accretion disc models, using the a-prescription, and local 
disc structure computations, the values of l/cr lie more in the range 
3-10 (FigureHl. 

Thus the simple theoretical disc models appear to display disc 
ages which are too high by factors of order 3-10. This remains 
true, even when late or delayed infall onto the discs is taken into 
account. 

Given this, we have investigated three further physical effects 
which can have a bearing on the measured disc age. These are 

(i) It may be that those parts of the disc which are sufficiently 
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Figure 9. Distributions of stellar age, accretion rate and disc mass of the observations compared between those systems which fall above and below the line 
idisc = The corresponding KS-test P- values are 0.954, 0.049 and 0.003 respectively. Both the stellar age and accretion rate distributions seem similar for 
high and low disc ages, but the are significantly more low mass systems at low disc ages, and more high mass systems at high disc ages. This may indicate 
there is a systematic error involved in the disc mass measurements. 



cool that the MRI cannot operate, and so cannot provide a viscos- 
ity. In this case, accretion initiall y occurs as a series of outbursts 
jArmitage, Livio & Pringlell200lh . This has a profound impact on 
the locus the disc follows in the — tdisc plane. The locus in 
this case is capable of covering the area in the plane occupied by 
observed systems (Figure |6}. However, in these models the discs 
spend most of their time at very large values of tdisc, with only 
brief sorties to much lower values (during the outbursts, when the 
accretion rates are high). This model could therefore account for 
the observations, but there would need to be very strong selection 
effects against the detection of systems with low accretion rates. 

(ii) In the standard models all the disc material is assumed to 
end up accreting onto the central star. In reality, it is possible that 
at late stages the disc material is blown away by the central star. 
The effect of this is to reduce the accretion rate onto the star, at late 
times, and so to correspondingly increase the disc age still further. 

(iii) One of the expectations of these 'proto-planetary' discs is 
that at least some of them succeed in forming planets. Formation of 
a planet in the disc which is large enough to interrupt disc accretion, 
and to allow migration of the planetary orbit, also has the effect of 
decreasing the central accretion rate, and so of increasing the disc 
age. 

It is therefore difficult for current models of the evolution of proto- 
planetary discs to provide a satisfactory fit to the observed points, 
taken at face value, in the t, — fdisc plane. Hence it is worth en- 
quiring whether some systematic effect or effects might lead one 
or more of the observational quantities to be mis-estimated. The 
three relevant quantities are (i) stellar age, (ii) disc mass, and (iii) 
accretion rate. 

If one of these quantities was primarily responsible for the 
observed disc ages to appear artificially low, it might be that its 
estimated value varies systematically with disc age. Thus we have 
investigated separately how each of the distributions of stellar age, 
accretion rate and disc mass compare, between those data points 
that lie above and below the tdisc = line. From Figure|9]it is ev- 
ident the distributions of stellar age and of accretion rate for those 
objects with high and low ratios of tdisc/t* are not significantly 
different. This is confirmed by the 2-sided K-S test on the distribu- 
tions for high and low disc age relative to stellar age (points either 
above or below the line), which returns probabilities that these are 
drawn from the same distribution of 0.954 and 0.049, respectively. 
However, the disc mass distribution does seem to show significant 
dependence on the value of fdisc/i*, and the K-S test indicates 
the probability of only 0.003 that the disc masses of objects with 



tdisc > i* and fdisc < i* are drawn from the same distribution 
(shown graphically in Figure [9](. Of course because disc age is an 
estimator of stellar age and is dependent on both Afdisc (linearly) 
and M (inversely) all of these quantities are interdependent. This 
clouds the meaning of the plots in Figure |9] somewhat, but what is 
clear is that the value of A/disc is dependent on whether the system 
has a high or low disc age relative to its stellar age, where as the 
values of M and are not. In particular, objects with low values 
of tdisc/t* tend to have on average lower values of Afdisc- 

To explain the discrepancy through either the stellar age or 
the accretion rate, one or both of these quantities would need to be 
systematically overestimated. Although the K-S tests hint towards 
the disc masses being the culprit of the discrepancy we will discuss 
the potential for mis-estimation of the other variables in the disc 
age. 

Measurement of accretio n rates contains several potential 
sources of systematic e rrors jCalvet, Hartmarm & Stron] l2000l : 
ICalvet & D'Alessiol201ll) . To measure the UV excess the emission 
from the star must be deducted, but stellar models do not always fit 
well to young stars. A bolometric correction must be applied to cal- 
culate the accretion luminosity, the error here could be up to a factor 
of 2 if line emission dominates over the continuum. The conver- 
sion between flux and luminosity also requires the distance to the 
star to be known. Finally, converting to the accretion rate requires 
knowledge of the mass and radius of the star and an assumption re- 
garding inner disc geometry. These errors are both significant and 
inter-dependent, and it is not clear if their combined effect will be 
in any particular direction or have an overall random effect on a 
large sample. 

Ages of young stars are measured by fitting theoretical pre- 
main-sequence evolutionary tracks to data in the H-R diagram, the 
position of which is derived by determining the stellar luminos- 
ity and its spectral type. The process is clouded by the presence 
of an accretion disc which introduces complications in determin- 
ing the correct luminosity (see ^2.\\ . this uncertainty is then in- 
creased as there is no consensus as to which theoretical track is the 
best fit, hence different au thors derive differen t age estimates, even 
for the same initial data. [Hillenbrand ('2009) concluded that this 
met hod likely und er-predicts low mass stellar ages by 30 ~ 100% 
and lNavloJ ( |2009|) recently developed a new fitting procedure, de- 
riving stellar ages of young clusters that were a factor of 1.5-2 
longer than those commonly found, suggesting that current age es- 
timates may be systematically too short. Thus, the errors in stellar 
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ages are likely in the opposite direction to that which would assist 
in explaining the discrepancy we have found. 

We therefore speculate that there might be a systematic error 
in the measurement of disc masses. We note that if all disc masses 
were currently being systematically underestimated by a factor of 
3-5, then we would achieve approximate agreement between the 
predictions of theoretical models and the observations. We show 
this in Figure [Tol 

The suggestion that there is a tendency to u nderestimate the 
masses of proto-planetary discs is not a new one, iHartmann et alj 
( l2006l) suggests that grain growth and the resulting reduction in 
dust opacity could lead to a systematic underestimation of disc 
masses. The basic problem is that disc mass is estimated from mea- 
surement of far-IR emission from dust grains. There may be a ten- 
dency for the gas to be evaporated, or otherwise lost, from discs 
at a greater rate than the grains. This would lead to disc masses 
being overestimated. Another problem is that measurement of disc 
mass relies on assumptions about the relative number of smaller 
grains (which dominate the radiative emission) and the larger grains 
(which dominate the grain mass). If grains tend to coagulate with 
time (a process which is a necessary precursor to the formation 
of planets) then the ratio of larger to smaller grains grows with 
time, a nd disc masses would be underestimated. I Williams & Ciezal 

If) suggest that the latter dominates, and that the overall ten- 
dency is that dust masses are underestimated. 
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